Astronomy & Astrophysics manuscript no. paper 


© ESO 2009 


June 24, 2009 





The onset of star formation in primordial lialoes 

Umberto Maio^'^, Benedetta Ciardi\ Naoki Yoshida^, Klaus Dolag' and Luca Tornatore^ 



o 
o 

(N 
(N 

u 

O 

(N 
> 
oo 

o 

\6 
o 
o\ 
o 

> 



' Max-Planck-Institut fiir Astrophysik, Karl-Schwarzschild-StraBe 1, 85748 Garching bei Miinchen, Germany 

^ Max-Planck-Institut fiir extraterrestrische Physik, GiessenbachstraBe 1, 85748 Garching bei Miinchen, Germany 

^ IPMU, U-Tokyo 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8568, Japan 

* INAF - Osservatorio astronomico di Trieste, Via Tiepolo 11, 34143 Trieste, Italy. 

Received ... 

ABSTRACT 

Context. Star formation remains an unsolved problem in astrophysics. Numerical studies of large-scale structure simulations cannot 
resolve the process and their approach usually assumes that only gas denser than a typical threshold can host and form stars. 
Aims. We investigate the onset of cosmological star formation and compare several very-high-resolution, three-dimensional, N- 
body/SPH simulations that include non-equilibrium, atomic and molecular chemistry, star formation prescriptions, and feedback 
effects. 

Methods. We study how primordial star formation depends on gas density threshold, cosmological parameters, and initial set-ups. 
Results. For mean-density initial conditions, we find that standard low-density star-formation threshold (0.2 h- cm"'') models predict 
the onset of star formation atz~ 25-31, depending on the adopted cosmology. In these models, stars are formed automatically when 
the gas density increases above the adopted threshold, regardless of the time between the moment when the threshold is reached and 
the effective runaway collapse. While this is a reasonable approximation at low redshift, at high redshift this time interval represents 
a significant fraction of the Hubble time and thus this assumption can induce large artificial offsets to the onset of star formation. 
Choosing higher density thresholds (135 cm"') allows the entire cooling process to be followed, and the onset of star formation is 
then estimated to be at redshift z ~ 12 - 16. When isolated, rare, high-density peaks are considered, the chemical evolution is much 
faster and the first star formation episodes occur at z > 40, almost regardless of the choice of the density threshold. 
Conclusions. These results could have implications for the formation redshift of the first cosmological objects, as inferred from direct 
numerical simulations of mean-density environments and studies of the reionization history of the universe. 

Key words. Cosmology: theory - early structure formation 



1. Introduction 

Understanding primordial structure formation is one of the fun- 
damental issues of modern astrophysics and cosmology. There 
is wide agreement that not only consists the universe of ordi- 
nary "baryonic" matter but also a large fraction of unknown 
"dark" matter, whose effects are only gravitational. Baryonic 
matter appears to constitute only a small fraction of the total 
cosmological matter content with a present-day density param- 
eter Qp f, = 0.0441 compared to Qo,m = 0.258 (Hinshaw et al. 
1200 8). Since the universe is observed to have zero curvature, i.e. 
to have a total density parameter Q.ojoi - 1 , these data imply that 
an additional density term exists Qq.a - 0.742. This is probabl y 
related to the so-called "cosmological constant" ('Einstein'1917'), 
or, as i nitially suggest ed by Ratra & Peebles ( 1988), Wetterich 
( Il988h . iBrax &Martid (119991) and IPeebles & Ratral (l2003l) . to 
other kinds of unknown "dark energies", whose effects on early 
struc ture formation history have been stud ied by e.g. Maio et al. 
(1200 6) with numerical simulations and bv lCrociani et al.l (I2008D 
with analytical calculations. 

The existence of non-baryonic matter was suggested several 
decades ago, and structure formation models based on the 
growth of primordia l gravitational instabilities ( lPeebleslll974l: 
IWhi te & Rees 1978) were developed following the early work 
bv lOunn & Go tt (1972). 

Hydrodynamical simulation codes [the first dating back to 



lEvrardl (Il988l) and iHernquist & Kat3 (Il989l) l have become 
a powerful tool, but because of computational limitations, 
plausible subgrid models have always been r equired to take 
into account star formation events (e.g. Cen & Ostrikei 



T9 92t lKatzlll992t iKatz et al.lfT996l: ISoringel & Hernauist 20031 
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Dal la Vecchia & Schavd l2008t ISchave & Dalla Vecchial l2008l) 
These simulations model the converging gas infall into dark- 
matter potential wells, by following the gas that becomes shock 
heated and subsequently cools by atomic and/or molecular cool- 
ing. Given the many orders of magnitude (in scale and density) 
spanned, it is computationally extremely challenging to simulate 
the process down to the formation of single stars. 
The gas physics in structure formation simulations has been 
typically approached with either lagrangian smoothed particle 
hydro-dynamics (SPH) or Eulerian mesh codes. A particular 
subclass is constituted by adaptive mesh refinement (AMR) 
codes, which allow further decomposition of the mesh around 
high-density regions, achieving a higher resolution. The main 
advantage of the SPH approach is its ability to follow self grav- 
ity in detail, while hydrodynamical instabilities are usually cap- 
tured by mesh codes. 

To account for star-formation episodes, both SPH and mesh 
schemes rely on spe cific assumptions. The prescriptions applied 
in mesh codes (e.g. ICen & OstrikeJ I 19921: llnutsuka& Mivamal 
IT99I iTVuelove et al.lll997l)" usually assume that the star forma- 
tion rate is proportional to th e density of overdense gas, while 
those used in SPH codes (e g. lKatzlll992l:lBate & Burker^l 19971: 
ISpringel & Hernquistll2003l) are based on the existence of a den- 
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sity threshold above which the gas is gradually converted into 
stars. Here we make use of SPH simulations. 
The typical timescales involved in the process of gas condensa- 
tion are the free-fall time, tjf, and the cooling time, tcooi- Gas 
condensation is expected to take place only if tcooi < t/f- 
The free-fall time is defined as 



371 



(1) 



where G is the universal gravitational constant and p the den- 
sity of the medium; the numeric factor (37r/32)'^^ is exact for 
spherical symmetry only. The cooling time is defined as 



3 nkBT 



where n is the number density of the gas, kg the Boltzmann con- 
stant, T the temperature, and J^{T, rit) the cooling function (en- 
ergy emitted per unit time and volume), which is dependent on 
both temperature and number densities, n,, of the species con- 
stituting the gas. In the low-density limifl for two-body interac- 
tions, between particles x and y, £, can be written as 



£.{T, n f, Hy) - A(T)n^ny 



(3) 



with the quantum-mechanical function A(T) depending on the 
temperature of the species considered, n i and (see for exam- 
ple MaioetaL 2007). At r > 10"* K, the cooling is dominated 
by collisions of hydrogen atoms, which is the most abundant 
species in nature - about 93% in number fraction - and X. scales 
approximatively as «^ (we indicate with rifj the hydrogen num- 
ber density). 

The physical conditions in which the first structures form are 
characterized by a primordial chemical composition: mostly hy- 
drogen, deuterium, helium, and some simple molecules, e.g. H2 
and HD. 

The primordial sites in which the first stars form are thought 
to be small dark-matter haloes with masses ~ lO^M© - as ex- 
pected from predictions based on self-sim ilar gravitational con- 
densation and chemical evolution (e.g. iTegmark et al.l [19971 
iTrenti & Stiavellill2009) - and virial temperatures r,,,> ~ lO'* K. 
Once they are born, they illuminate the universe and mark the 
end of the "dark ages". The radiation propagates in the vicin- 
ity of the individual sources and the im pact on the subsequent 
structure formation dRicotti et al.l l2002allbl 1200 8) can be very 
significant leaving imp rints by a means of feedback effects (see 



significant leaving imp rmts by a meat 
ICiardi & Ferrarall2b05L for a review). 



The low virialization temperatures of primordial haloes are 
enough neither to excite nor to ionize hydrogen and the lack 
of any metals means that the gas can c ool and eventually forin 
objects only via molecular transitions ( Saslaw & Zipo vlll967t 
Peebles & Dic"ka[T968t iHollenbach & McKed ll979: Maio et alj 



20071 for a detailed study of the cooling efficiency in differ- 



ent regimes). They have rotational energy separations with ex- 
citation temperatures below 10"* K, and therefore it is possi- 
ble to collisionally populate their higher levels with the conse- 
quent emission of radiation and resulting gas cooling. Since the 



' This widely-used approximation is appropriate as, according to the 
classical spherical "top-hat" model, a virialized object has a total mass 
density of l&n^ times the critical density, which corresponds, on aver- 
age, to a total number density of ~ 2h^ cm"' at z ~ 15, for a WMAP5 
cosmology and a mean molecular weight yu ^ 1. The transition to a 
high-density statistical equilibrium regime happens at critical number 
densities of ~ 10"* cm"'. 



molecular energy-state separations are typically smaller than the 
atomic ones, cooling will of course be slower, but still capable of 
bringing the temperature down to ~ 10^ K ( Yoshida et al.ll2003l; 



lOmukai & PallaH2003tr^hida et al.lll006tlGao et al.ll2'007ir 
To follow the entire process of structure and star formation 
in numerical simulations, one should implement the entire set 
of chemical reactions and hydrodynamical equations and from 
those calculate the abundance evolution and the corresponding 
cooling terms. In practice, performing these computations is 
very expensive and time consuming and it becomes extremely 
challenging to follow the formation of structures from the initial 
gas infall into the dark-matter potential wells to the final birth of 
star s. Neverthele s s, efforts are being made in this direction (e.g . 
A bel et al.l l2002t iBromm & Larsod l2004t lYoshida et al.] 120071; 



(2) iWhalen et al.ll2008h 



For this reason, more practical, even if sometimes coarse, sim- 
ple models are adopted. In brief, star formation relies on semi- 
empirical and numerical recipes based on chosen criteria to 
convert gas into stars and obtain the star formation rate, care- 
fully normalized to fit observational data at the present day. 
In particular, in SPH approaches a single particle represents a 
population of stars with assigned mass distribution. The stan- 
dard method used is to assume that once the gas has reached 
a given density thres hold it automat ically for ms starfl [e.g 
Cen & Ostriker ( 1992), Katz ( 1992), Katz et al. (1996) and the 
popular Springel & Hernquist ( 2003) model, inspired by the pre- 
vious works], regardless of the time between the moment when 
the threshold is reached and the effective run-away collapse, 
which typically takes place at densities ~ 10^-10"* cm"-^. While 
this is a reasonable approximation at low redshift, in "average" 
regions of the universe at high redshift this time interval repre- 
sents a significant fraction of the Hubble time and thus the as- 
sumption can induce large artificial offsets on the onset of star 
formation and influence the evolution in the derived star for- 
mation rate. Thus, extrapolations to high redshifts of the low- 
density thresholds (few 1"^ cm"^) used to model the star forma- 
tion rate in the low-redshift universe, may not always be justi- 
fiable. For this reason, high-redshift applications require higher 
resolutions and a higher density threshold. 
In this paper, we are interested in modeling star formation as a 
global process in regions of mean density in the universe, not 
directly in the very first stars, which instead form in highly over- 
dense, isolated regions. In particular, we discuss the importance 
of the choice of the density threshold for star formation in sim- 
ulations of early structure formation. We present a criterion to 
choose this threshold (Section |2]i and some test cases based on 
high-resolution simulations (Section [3] and HJ. Then we present 
our results and conclusions (Section|5]l. 

2. Threshold for star formation 

Accor ding to the u sual scenario of structure formation, the Jeans 
mass (|Jeanslll902l) is the fundamental quantity that allows us to 
distinguish collapsing from non-collapsing objects, under gravi- 
tational instability. For a perfect, isothermal gas, it is given by 



M 



6 [finiHGl 



P 



-1/2 



(4) 



where mn is the mass of the hydrogen atom, and T and p are 
the temperature and density of the gas, respectively. At very 



- To reduce the computation time for calculations of fragmentation 
sometimes, part icles with densit ies above the threshold are replaced by 
'sink' particles jBate et aljl995t) . 
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high redshift (z ~ 30 - 20), typical haloes have masses of 
~ 10^ - 10*" Mq, which can increase up to ~ 10^ - 10'^ Mq by 
z~ 10. 

As mentioned in the introduction, the density threshold for 
star formation in numerical simulations is typically fixed to 
some constant value, irrespective of the simulation resolution. 
However, it would be desirable to have a star formation criterion 
that allows us to reach scales that fully resolve the Jeans mass. 
The SPH algorithm implicitly imposes a minimum mass res- 
olution limit, because to compute the different physical quan- 
tities, a fixed number of neighbours (i.e. number of particles 
within the smoothing length|3 , h) is used. This also induces a 
minimum resolvable mass, which is the total mass of neigh- 
bouring particles. This is particularly important in SPH sim- 
ulations of cosmological structures and galaxy formation, be- 
cause only if the minimum resolvable mass is far smaller than 
the Jeans mass, it is possible to ensure that the results are not 
affected by numerics nor by th e details of the implementation 
adopted dBate & Burkerl|[l997h . Otherwise, unresolved, Jeans 
unstable clumps can easily be found to exhibit unphysical be- 
haviour (e.g. over- fragmentation problem in l ow-resolution sim- 
ulations). Furthermor e, it was shown (Navarr o & Whitdll993l; 
iBate & BurkertI [T997I) that the minimum number of particles 
needed to obtain reasonable and converging results is about 
twice the number of neighbours (~ 10^ particles). 
If Mres is the gas mass resolution of a given simulation, we can 
assume that: 



Mj = NMres, 

where » 1 and that the critical threshold is 

7 / ^ n3 



Pth - 



knT 



1.31 ■ 10-1'' IM,,, 



(5) 

(6) 
(7) 



For Mres = 10^ Mg, T = 10^ K, and ju = 1 and when 
using A' = 10^ gas particles, one has p,/, ~ lO^^'gcm"-', 
corresponding to a physical number density of ~ 10^ cm"^. 
The above equations should be considered as a guideline 
to estimating the density t hr eshold . We note that we largely 
fulfill the .Bate & Burkerti ( Il997h conditions, because the 
neighbour number in our simulations is Nneigh - 32 (see 
Section |3), so we resolve the Jeans mass with about three 
times the number of neighbours (A^ ^ 3N„eigh)- As akeady 
mentioned however, commonly adopted density thresholds 
are for practical reasons usually chos e n to be of order of 
-10- 
200l 



cm ^ or less ( Katz et ^ 1996t Springel & HernauistI 



IScannapieco et al. 12005 : 'Governato et al.' 120071: 

Tornato re et al. '2007b; Schave & Dalla Vecchia 20081). For 



of 10" cm" , while the gas-particle mass was close to be the 
Jeans mass and the gravitational softening of order of the Jeans 
length. 

In general, in simulations of both cosmic structure and galaxy 
formation, it is quite hard to fully resolve the Jeans mass of 
the collapsing fragments, and star formation is often assumed 
to occur while the gas fal ling into the dark-ma tter potential 
wells is still heating up. The lBate & Burker3 (1 19971) requirement 



^ The definition of smoothing length may vary from authors to au- 
thors. It is sometimes meant to be the width of the SPH smoothing ker- 
nel, but other times the length scale on which the SPH smoothing kernel 
becomes zero. 



is not usually satisfied, since the main goal is usually not to 
follow the entire process of collapse and fragmentation, but to 
obtain a qualitatively representative sample of the cosmological 
evolution. 

A high value for the threshold is also important to capture 
the relevant phases of cooling. In the following, we show that 
molecule radiative losses, at temperatures of ~ lO-' - 10^ K 
where they can balance the heating of the infalling gas, produce 
an isothermal state and a subsequent cooling regime. For a re- 
gion of mean density, the time spent by the gas in the isothermal 
state can be a substantial fraction of the Hubble time. Therefore, 
it is important for the threshold to be on the right-hand side 
of the peak in the phase-diagram (i.e. at densities higher than 
the isothermal regime), so that the delay between reaching the 
threshold and the true star formation is negligible (see Section 
4). 

In the following, we investigate the effect of different choices 
of star formation thresholds at high redshift, describe the 
simulations performed, and discuss the results obtained. 

3. Simulation set-up 

To study the effect of different threshold prescriptions on the 
onset of star formation, we completed very high resolution, 
three-dimensional, hydrodynamic simulations including non- 
equilibrium atomic and molecular chemistry, star formation, and 

wind feedback. 

We used the code Gadget-2 (ISpringelll2005h in its modified 
form, which includes stellar evolution and metal pollution 
(iTornatore et al.l l2007ah . primordial molecular chemistry (fol- 
lowing the evolution of e", H, H+, He, He+, He++, H2, H+, H", 
D, D^, HD, HeH^), and fine structure metal transiti on cooling 
(O, C\ Si\ Fe^) at temperatures lower than 10'* K (Maio et aJJ 
I2OO7I I2OO8I). We perform hydro-calculations by fixing the num- 
ber of SPH neighbours to Nneigh = 32. 

The simulations have a comoving box size of L = 1 Mpc and 
sample the cosmological medium with a uniform realization of 
320^ particles for both gas and dark-matter species (for a total 
number of 2 x 320^ particles). The resulting gas-particle mass is 
of the order of 10^ Mq, which is consistent with the discussion 
in the previous section. 

We note that this configuration enables us to easily resolve the 
Jeans length for shock heated/cooling cosmic gas (the Jeans 
length for gas with T ~ 10* K and p ~ 10"^^ - 10"^'*g/cm-' is 
~ 3 - 1 kpc, much longer than the comoving gravitational soft- 
ening of ~ 0.1 kpc). 

The initial conditions (set at redshift z - 100) are generated with 
a fast Fourier transform grid of N,„esh = 320 meshes and a max- 
imum wave-number (Nyquist frequency) 



example, Wi ersma et alJ (l2009l) used a number density threshold ^ 



InN, 



Nyquist 



mesh 



2L 



1 kpc 



(8) 



(i.e. a minimum wavelength of ILjNmesh - 6.25 kpc), so that, 

for each wave-number, ||k|| < kf^yqmst- 

We will refer to this sampling as "mean region". 

We considered two different sets of cosmological parameters: 

- standard model: Qo,m = 0.3, Qo,a = 0.7, Qo,* = 0.04, h — 
0.7, cTg = 0.9 and n - I, where the symbols have the usual 
meanings. The coiTesponding dark-matter and gas-particle 
masses are ~ 755 Mq/Zi and ~ 1 16 M0//1, respectively. 



The comoving gravitational softening is usually estimated as 1 /20 
or 1/30 the mean inter-particle separation. 



4 



Umberto Maio et al.: The onset of star formation 



Table 1. Parameters adopted for the simulations. 



ivioQei 


Number of 












u 


C 8 


n 


SF threshold 




gas+dm particles 


[Mo//!] 


[Mo//!] 
















wmap5-ht 


2 X 32768000 


128 


621 


0.258 


0.742 


0.0441 


0.72 


0.8 


0.96 


135.0 


wmap5-lt 


2 X 32768000 


128 


621 


0.258 


0.742 


0.0441 


0.72 


0.8 


0.96 


0.2 


std-ht 


2 X 32768000 


116 


755 


0.300 


0.700 


0.0400 


0.70 


0.9 


1.00 


135.0 


std-lt 


2 X 32768000 


116 


755 


0.300 


0.700 


0.0400 


0.70 


0.9 


1.00 


0.2 


zoom-std-ht 


2x41226712 


3.9 


25.6 


0.300 


0.700 


0.0400 


0.70 


0.9 


1.00 


135.0 



The columns (from left to right) specify: name of the run, number of particles used, gas-particle mass, dark-matter-particle mass, Qo.m, ^o a, i^o.*, 
h, (Tg, spectral index, star formation density threshold. 



- WM AP5 model: da ta from 5-year WMAP (WMAP5) satel- 
lite (iHinshaw et alJ '2008) suggest that Qo,,, - 0.258, Qqa = 
0.742, Qo,z> = 0.0441, h = 0.72, erg = 0.796, and « = 0.96. 
In this case, the coiTesponding dark-matter and gas-particle 
masses are ~ 621 M^/h and ~ 128 Mq/Zi, respectively. 

Following the discussion in the previous sections, we also con- 
sider two different values for the star formation density thresh- 
old: 



a low-density threshold of 0.2/z^cm"^ (physical), compat- 
ible with the one adopted in the Gadget cod e and the 
ones widely used in the literature ( for example Katz et alj 
1996t ISpringel & Herng uist 2003t IXornatore et all 12007^ 



Pawlik et al 



20091) : 



- a high-density threshold of 135/!^cm"^ (physical), as com- 
puted from Eqs. (|6]l and (Q. This value is adequate for mod- 
elling atomic processes even in small ~ 10^ Mq haloes at 
z ~ 20. Moreover, this threshold typically falls in density 
regimes where cooling dominates over heating, allowing us 
to properly resolve gas condensation down to the bottom of 
the cooling branch. 

A summary of all the simulation features is given in Table[T] We 
denote with the labels "std" and "wmap5" the runs with standard 
and WMAP5 cosmology, respectively, and with "It" and "ht" the 
runs with low- and hig h-density thresholds, re spectively. 
We note that the Springel & Hernqiiisj (|2003|) model used here 
to describe the star formation process is strictly applicable only 
as long as more than one star per SPH particle is present, i.e. each 
SPH particle is considered as a 'simple stellar po pulation' with 
a given mass distribution. Although some studies ([Yoshida et all 
120031: iBromm et al.l l2002t iBromm & Larson 2004) seem to in- 
dicate (or assume) that the very first episode of star forma- 
tion could result in a single, very massive star per halo, this 
should apply preferentially to very high redshift, high density, 
isolated objects. In any case, the exact shape of the IMF of 
primordial stars ( e.g. Schwarzschild & Spitzer 1953; Larson] 
Il998t iNakamura & Umemura .2001.: Omukai & Palla 2003;) is 
still a matter of speculation and lively debate. For this reason 
and because we are interested mainly in the global star forma- 
tioi i process, we ii s ed the [Springel & Hernguist (2003) model 
[as iTornatore et al.l (l2007bl) also did to describe both a primor- 
dial top-heavy and a more standard star formation mode] and al- 
lowed the IMF to be a free parameter (although in the test cases 
reported here we always adopt a top-heavy IMF). The aim of this 
paper is to investigate the effects of the density threshold on star 
formation, and we leave discussion on the IMF to future work. 
Finally, to investigate primordial star formation events in local 



high-density regions, we perform a very high-resolution numer- 
ical simulation of a rare high-sigma peak with comoving radius 
~ 140kpc//!. This region is selected using the zoomed initial 
condition techniqu e on a ~ 10^ M p halo formed in a dark-matter- 
only simulation ( Gao et al.l2007lF l We divided each particle into 
gas and dark-matter component, according to the standard model 
parameters. The resulting gas-particle mass is ~ AM^/h and 
dark matter particles have a mass of ~ 26 M^/h (in Table [1] this 
simulation is labelled "zoom-std-ht"). 

By a quick comparison of the different parameters adopted, we 
expect that, once the density threshold has been fixed, the stan- 
dard cosmological mean-region simulations will show earlier 
structure formation episodes with respect to the corresponding 
wmap5 ones. This is because they have higher spectral parame- 
ters and higher matter content. The high-density region is a bi- 
ased over-dense region already at early times, and therefore, its 
evolution is expected to be much faster 

4. Results 

We present the results of the simulations with the sets of param- 
eters described above. We discuss first the mean region of the 
universe (Section |4Tt and then the high-density region (Section 



4.1. Mean-region simulation 

Our reference run is the wmap5-ht model with init i al co rtipo- 
sition given by the values quoted in iGalli & Pallal ( Il99 8n at 
z - 100. We show some evolutionary stages in Fig. [1] (upper 
set of panels). In the maps, the first column refers to tempera- 
ture, the second to gas density and the third to molecular frac- 
tion at z = 30.16 and z - 12.17, respectively. The creation 
of new molecules is clearly evident, together with the related 
growth of structures. More specifically, as time passes, one can 
see the heating undergone by the gas in dense regions, because 
of structure formation shocks. The temperature increases from a 
few hundreds Kelvin in the low-density regions, to ~ lO'^K in 
the denser regions. In the meantime, the molecular fraction also 
evolves accordingly up to values higher than lO"**. Soon after, 
the production of molecules increases rapidly (up to ~ 10"^) aid- 
ing the star formation process, which, for this simulation, starts 
atz ^ 12. 



' We use the "R4" initial conditions presented there. 

* We assume a primordial neutral gas with residual electron and 
fractions x^- ^ xh+ - 4 • 10""*, H2 fraction xh, = 10"*, H2 fraction 
xh+ = 3 ■ 10-2', D fraction xq = 3.5 • lO-\ HD fraction xhd = 7 • 10-'^ 
D+ fraction xd+ = 4 ■ 10"', HeH+ fraction XHeH+ = lO"'". 
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Fig. 1. First, second, and third column are respectively temperature, density, and molecule maps. The first two rows refer to the mean-region 
simulation at redshift 12.17 (top) and 30.16 (bottom). The box size is 1 Mpc comoving. The last two rows refer to the high-density region at 
redshift 50 (top) and 70 (bottom). The region size is ~ 140 kpc/Zi comoving. All quantities are smoothed on a 276-pixel side grid. 
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Fig. 2. Upper panel: phase diagram at redshift z = 12.17 (just before 
the onset of star formation) for the wmap5-ht simulation. The verti- 
cal straight lines indicate a low physical critical density threshold of 
0.2 A^cm"^ (dashed line) and a higher physical critical density threshold 
of 135 /j-cm"' (solid line). Lower panel: average effective index com- 
puted over the whole range of densities. The three horizontal dotted 
lines show values of 5/3, 1 and 1/3, respectively from top to bottom. 
The solid line shows a and the dashed line shows y (see text for defini- 
tions). 




10^ r 



10 



1 2 



Fig. 3. Evolution of the ten most massive haloes in the wmap5-ht cos- 
mological simulation (dotted lines). The redshift at which the first star 
forms in each halo is indicated by the filled star symbols. After that, star 
formation continues along the solid lines. 



In Fig.|2] we show the phase diagram (comoving density versus 
temperature) at redshift z ^ 12, i.e. just before the onset of star 
formation. The low-density gas, which is shock-heated by the 
collapse of the first primordial haloes, is seen on the left side of 
the panel. Starting from values for a temperature of ~ 10^ K, the 
gas is progressively heated to ~ lO'* K and moves along the ris- 
ing branch. At this stage, collisions become more frequent due 
to the higher temperature. The upper energy levels of particles 
become excited and the subsequent de-excitation is accompa- 
nied by the emission of radiation. This effect is negligible at low 
densities, because collisions are rare and the fraction of energy 
converted into radiation is small. When the density increases, 
the cooling becomes comparable to the heating and an isother- 
mal regime with no significant net change in the temperature 
is reached. This appears at the tip of the phase diagram (and 
in the behaviour of the effective index, as discussed below), at 
r ~ 10"* K, where the cooling is dominated by atomic Lye tran- 
sitions and accompanied by runaway collapse. At higher densi- 
ties, radiative losses overtake heating and induce a fast cooling 
phase (dominated by molecules, mostly H2). 
The solid vertical line corresponds to the physical high-density 
star formation threshold (135 /j^cm"^) and, for comparison, we 
also plot the dashed line for a physical number density of 
0.2/z^cm"^. We stress that by adopting a low-density thresh- 
old for star formation one completely misses the isothermal and 
cooling part of the phase diagram, and thus a correct modeling 
of the cooling regions within the simulations. This can affect the 
onset of star formation, particularly at high redshift, when the 
time needed for the gas to evolve from the low-density threshold 
to the high-density threshold (~ 2 • 10^ yr) can be a substantial 
fraction of the Hubble time (~ 4 x 10^ yr at z ~ 12). We note 
that the time elapsed between the attainment of the isothermal 
peak in the phase diagram and the end of the cooling branch is 
~ 6 X 10' yr. The evolution that follows the end of the cooling 
branch is characterized by the format ion of a dense core, which 
accretes gas on free-fall timescales ( Yoshida et al.l i2006). This 
phase has a very short duration (~ 10^ - lO^yr) during which 
the central densities increase to ~ lO'^cm The problem is 
less severe at lower redshift, when the Hubble time becomes of 
the order of several Gyr 

The density and temperature behaviour can also be described by 
an effective inde?(Q, which depends on the physical conditions of 
the gas regime considered. In the lower panel of Fig.|2] we plot 
the effective index as a function of density. The solid line refers 
to the value a = \ + {dT IT)j(dplp), which takes into account 
changes in the sign of the temperature derivative, distinguishing 
the heating regime {a > 1) from the cooling regime (a < 1) . 
The dashed line refers to y = 1 + \{dT IT)l{dplp)\, so that y is al- 
ways > 1 . Dotted horizontal straight lines INDICATE values of 
5/3, 1, and 1/3. In correspondence with the isothermal peak in 
the T -p plane, it is o- = 7 = 1, which marks the transition from 
the heating to the cooling regime. At this stage, we expect the 
gas runaway collapse to begin and last for the following cooling 
regime, at which point a oscillates around the value of 1/3. 
In Fig. [3] we plot the evolution of the ten most massive haloes 
found in the simulation. We also show the redshift at which stars 
are produced (filled star symbols) in each object. The haloes are 
found using a friend-of-friend algorithm with a linking length 
equal to 20% of the mean inter-particle separation. Typical halo 
masses at redshift z ~ 12, when star formation starts , are of the 
order of W (see also IWise & Abelll2007l lIOOSi) and reach 



' By effective index of the gas, 7, we mean P cc , with P pressure 
and p density. This is simply related to the politropic index. 
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Fig. 4. Star formation rate as a function of redshift for the differ- 
ent models, from left to right: WMAP5 cosmology and high-density 
threshold (solid red line), standard cosmology and high-density thresh- 
old (dotted blue line), WMAP5 cosmology and low-density threshold 
(dot-dashed black line), standard cosmology and low-density threshold 
(long-dashed-short-dashed magenta line). The green short-dashed line 
refers to the simulation of the high-density region with standard param- 
eters and high-density threshold. 



densities of ~ 10^ cm"^. 

For comparison, we completed the same simulation using stan- 
dard cosmological parameters (std-ht run). In this case, the over- 
all picture is similar, but we detected a faster evolution, with 
earlier structure formation, as expected from the previous dis- 
cussion in Section[3] The first star formation events are detected 
at redshift z ~ 16 in haloes with masses ~ 10^ Mq. 
This can clearly be seen in Fig.|4] where we plot the star forma- 
tion rate as a function of redshift for the different simulations (to 
compute the star formation rate, we adop t the implementation 
described bv lSpringel & Hernquis3 (l2003l) ). 
The onset of star formation in the wmap5-ht model (red solid 
line) is delayed compared to the std-ht model (blue dotted line). 
For the wmap5-lt (black dashed line) and std-lt (magenta short- 
long-dashed line) models, star formation starts at z ~ 25 and 31, 
respectively. Thus, at these high redshifts, even small changes 
in the cosmology can be significant for the onset of star forma- 
tion. This is easily understood in terms of spectral parameters: 
the standard cosmology has higher spectral index and normaliza- 
tion; therefore, assigning more power on all scales with respect 
to WMAP5 values, leads to structure formation occurring much 
earlier. 

The choice of the density threshold makes an even larger dif- 
ference to the onset of star formation. In Fig.|4] the rates corre- 
sponding to the wmap5-lt (black dot-dashed line) and wmap5-ht 
(red solid line) show that star formation starts at z ~ 25 and 
12, respectively. The major difference between low- and high- 
density threshold models is that, in the former, the gas reaches 
the critical density much earlier. So, the redshift difference in the 
onset corresponds to the time that the gas needs to move from the 



low- to the high-density threshold (see Fig.|2]i. 
In addition, the simulations adopting the high-density thresholds 
slightly overtake the respective low-threshold cases. This hap- 
pens because the former did not remove the gas at higher red- 
shifts, it accumulated and ended in delayed bursts of star for- 
mation. Later, the star formation rates were restored to the same 
level. 

As already mentioned, the low-density threshold model is very 
commonly used both in numerical and semi-analytical works, 
because it does not require incorporation of molecular chemistry 
(the threshold being lower than the typical densities at which 
molecules become efficient coolants) and therefore it is easier to 
implement and allows faster simulations. However, it can com- 
promise the entire picture if the results are extrapolated to high 
redshift, when molecules are the main coolants and the time de- 
lay between the attainment of the low-density threshold and the 
bottom of the cooling branch occupies a significant fraction of 
the Hubble time. 

4.2. High-density region simulation 

We show results for the high-density region described in Section 
[3]and initialized at redshift z - 399. 

In this case, the physical number densities at the beginning 
of the simulation are in the range ~ 0.5 - SOh^cnT^ (at 
z ~ 200), with an average of ~ 4 cm"^, higher than the typical 
value adopted for the low-density threshold for star formation. 
Therefore, the conventional low-density model would produce 
unreasonable star formation at z ~ 200. To avoid this, it is 
common to add a further, additional, ad hoc constraint, which 
allows star formation only if the simulation over-densities are 
higher t han a given mini mum value - usually between ~ 50 and 
~ 100 ( Katz et al]|1996 , in Section 4.2, for example, suggest 
55.7). Thus, in this case it is this additional constraint that 
determines when the onset of star formation occurs, rather than 
the low-density threshold. 

We therefore ran a simulation with only a high-density thresh- 
old. For the sake of comparison, we still used the value of 
135/;^cm"^, although rigorously, following Eqs. (|6]l and (Q, 
one should adopt a value ~ 9 x 10^ cm"-' for a 3.9 M^/h gas- 
particle mass. Nonetheless, we checked that this choice does not 
affect our conclusions, since the threshold is already beyond the 
isothermal peak, in the fast cooling regime, where the timescales 
are extremely short (-10^ yr). All th e initial abundances are set 
according to the values suggested bv lGalli & Pallal (ll998lF l 
The simulation maps are shown in Fig. [T] (lower panels). As 
for the mean-density regions, they refer to temperature, density, 
and molecular fraction at redshifts z=50 and 70. As expected, 
we highlight that structure formation occurs far earlier than 
the mean-density case. Molecular abundances of ~ lO"*" are 
reached faster than for the mean-density region, where such a 
high fraction is found only at z < 30. Similarly, values of ~ 10"'* 
are already reached at z ~ 50 - 40, rather than z ~ 20 - see also 
discussion in Section|5]and Eq. (|9]l. 

In Fig.|5] we show the phase diagram and the behaviour of the 
effective index as a function of the comoving gas density at red- 
shift z - 45. Physical critical density thresholds of 0.2 h^cmT^ 
(dashed line), 135/i^cm"-' (solid line), and ~ 9 x h^croT^ 



The initial abundances (at z = 399) are consistent with a primordial 
neutral plasma having residual electron and fractions of Xe- - xh+ 
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Fig. 5. Upper panel: phase diagram at redshift z = 45.41 for the high- 
density region simulation for a purely non-equilibrium chemistry run 
(i.e. without star formation). The vertical straight lines show physical 
critical-density threshold of 0.2 /i^cm"' (dashed line), 135 li^ cm^^ (solid 
line) and 8.9 x 10"* h^cnr^ (dot-dashed line). Lower panel: average ef- 
fective index computed over the whole range of densities. The three 
horizontal dotted lines show values of 5/3, 1, and 1/3, respectively, from 
top to bottom. The solid line refers to a and the dashed line to y (see 
text for the definitions). 



(dot-dashed line) are marked in the figure. While the first two 
are the same as for the mean-density region simulation, the last 
one corresponds to the value obtained using Eqs. (|6]l and (|7]i. 
To emphasize the different characteristics of the phase diagram 
compared to the one obtained for the mean-density region, the 
plot was extended to densities higher than before. The isother- 
mal peak is reached at redshift z ~ 50. Unlike the mean-density 
simulation, the gas does not spend time on the isothermal 
plateau, but cools very rapidly (in less than 7 x 10^ yr) from 
~ lO-''^ K to ~ 10^ K and condenses into comoving densities of 
Pcom ~ 10"^' /i^g/cm-'. The rapidity of these events is reflected 
in the lack of particles in the intermediate stages of the cooling 
branch. 

As before, we also plot the effective gas index. The usual initial 
shock-heating behaviour and the following cooling is recovered 
up to much higher densities. At the bottom of the cooling 
branch, we find values of a that oscillate around 1/3 and 1. 
Since the last stages are quite fast, the low number of particles 
present introduces some statistical noise, which is evident in the 
plot. 

With our choice of the threshold, star formation sets in at z ~ 48 
(see Fig. |4]i. The additional time needed to reach the highest 
densities at the bottom of the cooling branch is extremely short 
(~ lO*" yr), so our choice ensures that the onset of star formation 
is correctly estimated. As there is no obvious, standard way of 
quantifying the star formation rate in these simulations, we do 
this by dividing the stellar mass formed at each time-step by the 
volume of the gas contained in the high-density region (a sphere 



of about 140 kpc//i radius). 



5. Discussion and conclusions 

We have studied the effect of different choices of the density 
threshold on the onse t of cosmic star formation in numerical 
SPH simulations (see iMaio et all 120071: iTornatore et al.ll2007al 
for technical details). 

In the literature, several studies are presented that follow the 
birth of primordial stars in early proto galaxies (examples are 
lYoshida etal.]|2006t IWise & Abelll2007h . These are mainly fo- 
cused on the initial phases of star formation and the effects on the 
immediate surroundings, which typically do not address more 
general issues such as the global star formation process (i.e. in a 
region at mean density rather than in high d ensity peaks), metal 
enrich ment, and IGM reionization (e.g. iBolton & Haehnela 
2007). Here instead, we are interested in simulating the more 
global star formation process in the high-redshift universe, cap- 
turing the relevant timescales and physical processes, i.e. the 
atomic and molecular physics that regulates the formation of pri- 
mordial stellar population. 

We have thus run simulations using initial conditions appropriate 
to a region of the universe with mean density and as a reference, 
using the zoom technique, a high-density peak. 
A basic process that leads to star formation, i.e. gas shock heat- 
ing up to ~ 10^ - 10^ K by infall into dark-matter haloes followed 
by radiative losses due mainly to molecular collisional excita- 
tions, is common to both scenarios. The main difference is asso- 
ciated with the global dynamics and timescales of the process. In 
the rare high-sigma peak, because of the higher densities, chem- 
ical reactions are faster and much more efficient with respect 
to the simulations of mean-density initial conditions. Therefore, 
the molecular fraction increases more rapidly, reaching a num- 
ber fraction of ~ 10"^ by z ~ 50 - 40 (compared to z ~ 20, 
for the corresponding mean-density case). These values are suf- 
ficient to make collisional cooling to dominate over heating and 
induce star formation episodes. 

Density and temperature behaviour can be described by an ef- 
fective index that depends on the physical conditions of the gas 
considered. Roughly, it is isothermal during the transition from 
the heating to the cooling regime, then it collapses (see the ef- 
fective index computed in the lower panel of Fig. |2]l until the 
bottom of the cooling branch is reached. More quantitatively, 
when cooling is dominated by H2, the cooling time in Eq. (|2]i 
can be approximated as 



kuT 



2 AH,(r) XH,nH 



(9) 



where x^^ is the H2 number fraction, AH,(r) is the H2 cooling 
function at temperature T, and the other symbols have their usual 
meanings. 

For gas at the beginning of the cooling branch, T ~ 10^'^ K and 
xiij ~ 10""*, giving tcooi ~ 7 X 10^ yr («h in cm"^). In the 
mean-density case (see phase diagram in Fig.|2]i, rin - 0.3 cm"-', 
while in the high-density region (see phase diagram in Fig. |5]l, 
hh - 6 cm"-'. This translates into a characteristic cooling time of 
~ 2 X 10^ yr for the former case and ~ 10^ yr for the latter 
These estimates show the relevance of following the full cool- 
ing branch when simulating star formation at high redshift in re- 
gions of mean density, because the characteristic cooling times 
are a substantial fraction of the Hubble time. This problem is 
less severe for simulations of high-density peaks, in which the 



Umberto Maio et al.: The onset of star formation 



9 



timescales are much shorter. 

If one considers a 10^ halo at z ~ 40, its average gas num- 
ber density is of the order of ~ 20//icm"^ (being /i the mean 
molecular weight), its virial temperature Tvii- ~ 7 x lO^juK and 
its cooling time tcoi ~ 3 x 10^ - 3 x 10^ yr for an H2 fraction 
of 10"^ - 10"^, respectively. The same halo at redshift z ~ 15 
would have an average gas number density of ~ l//icm"^, 
Tvii- ~ 5 X lO^K, and W ~ 7 x 10^ - 7 x 10'' yr. This means 
not only that objects of similar mass cool (and thus harbor star 
formation) on very different timescales according to their envi- 
ronment, but also that gas in high-redshift haloes can collapse 
and fragment much faster than in low-redshift ones, because its 
cooling capabilities are more efficient. 

For these reasons, though low-density thresholds can be 
useful tools to reproduce empirical surface-density rela- 
tions, such as the Kennicu tt-Sch midt law jK ennicutt 1998; 
ISpringel & HernquistI 120031: ISchave & Dalla Vecchial I2008D . 
physical insights can rely only on high-density thresholds. 
Indeed, imposing star formation events before the isothermal 
peak is reached could result in an artificially high-redshift for 
the onset of star formation. 

For the test cases presented in this paper, the value adopted 
for the high-density threshold is 135/i^cm"^, well beyond the 
isothermal peak of the gas. This allows a correct estimate of the 
relevant timescales, as the gas spends most of the time in the 
isothermal phase. In addition, following the evolution of the gas 
to higher densities allows highe r resolution of, e.g. , the morphol- 
ogy and disk galaxy structure dSaitoh et al.ll2008h . the dumpi- 
ness of the gas, and the features of the interstellar or intergalactic 
medium. On the other hand, running high-density threshold sim- 
ulations to the present age (z - 0) is computationally very chal- 
lenging because of the extremely short timescales involved in 
the calculations. Only simulations performed with a low-density 
threshold are currently run to z = and fine-tuned to reproduce 
the observed low-redshift evolution of the star formation den- 
sity. 

We note that the mean-density region, because of its small di- 
mensions, lacks massive haloes. In larger simulations, we expect 
to find rarer, larger, haloes, which can grow faster and host star 
formation in ~ 10^ Mq - lO'' Mq haloes. 

We have performed high-resolution, three-dimensional, N- 
body/SPH simulations including non-equilibrium atomic and 
molecular chemistry, star formation prescriptions, and feedback 
effects to investigate the onset of primordial star formation. We 
have studied how the primordial star formation rate changes ac- 
cording to different gas-density threshold, cosmological param- 
eters, and simulation set-ups. Our main findings are summarized 
in the following: 

- The typical low-density thresholds (below ~ 1 cm"^) are in- 
adequate for describing star formation episodes in mean re- 
gions of the universe at high redshift. To correctly estimate 
the onset of star formation, high-density thresholds are nec- 
essary. 

- In rare, high-density peaks, the density can be higher than the 
usual low-density thresholds from very early times, therefore 
these prescriptions are not physically meaningful. Density 
thresholds lying beyond the isothermal peak (several par- 
ticles per cm^, in our case ) are s till required: they should 
satisfy the iBate & BurkertI (1 19971) requirement (A^ at least 
~ 2N„eigh) and Eq. (|6]) of Section |2l However, as long as 
they are beyond the isothermal peak, given the faster evolu- 
tion in the phase diagram of the cooling particles in dense 
environments, the very exact value is not crucial. 



- Different values of the threshold and the cosmological pa- 
rameters can cause the onset of star formation at very differ- 
ent epochs: with a low-density threshold (0.2 cm"-'), star 
formation starts at z ~ 25 -3 1 (depending on the cosmology), 
while high-density threshold models (135 Z;^ cm"^) predict a 
much later onset, atz~ 12-16 (depending on the cosmol- 
ogy)- 

- Performing primordial, rare, high-density region simulations 
within the high-density threshold model, we find that the lo- 
cal star formation can set in as early as z ~ 48. 

We conclude by adding a few comments on the meaning of the 
term "onset of star formation". From a purely physical point of 
view, the onset of star formation occurs when the proton-proton 
nuclear reactions ignite in a collapsed, dense core. Nowadays, 
in numerical simulations of cosmic structure evolution, this def- 
inition cannot be adopted, since it is not feasible to follow and 
resolve the behaviour of the gas on the very large range of scales 
involved. Therefore, the onset of star formation is meant as the 
attainment of a certain density threshold: once SPH particles 
reach it, they are assumed to be dense enough to collapse and 
to host star formation. This threshold can be viewed as a point 
of no return, because it imposes a limit above which the natural 
gas evolution is strongly altered by star production and feed- 
back effects. The time when the threshold is reached and the 
time when the actual star formation takes place are typically not 
the same, but if the threshold is set appropriately, as discussed in 
the present work, they become very similar Broadly speaking, 
one could consider assigning a simple prescription to standard 
numerical simulations based on the typical delay time of gas in 
fall. In this way, the onset of star formation could be easily cor- 
rected without implementing high-density thresholds or molecu- 
lar evolution. In practice, this is not a trivial task: the time spent 
in the isothermal regime depends on numerous "environmental" 
factors, so different simulations with different initial conditions 
will be affected differently, according to their particular features. 
However, as an estimate, a typical free-fall time of ~ 10^ Myr is 
the one we expect in correspondence of ~ 10 ' cm"^ (conven- 
tional low-density thresholds). 

Throughout this paper, we have considered the onset of star for- 
mation to be the attainment of the density threshold, provided 
that the Jeans mass is resolved by a "large" number of SPH par- 
ticles and the isothermal peak in the phase diagram is resolved 
(and we have seen in Section |2] and |4] that the latter two condi- 
tions are strongly related). 
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